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ABSTRACT 

The  realizability  of  Reynolds  stress  models  in  homogeneous  turbulence  is  critically  as¬ 
sessed  from  a  theoretical  standpoint.  It  is  proven  that  a  well  known  second-order  closure 
formulated  by  Shih  and  Lumley  using  the  strong  realizability  constraints  of  Schumann  is,  in 
fact,  not  a  realizable  model.  The  problem  arises  from  the  failure  to  properly  satisfy  the  nec¬ 
essary  positive  second  time  derivative  constraint  when  a  principal  Reynolds  stress  vanishes 
-  a  fatal  flaw  that  becomes  apparent  when  the  non-analytic  terms  in  their  model  are  made 
single- valued  as  required  on  physical  grounds.  It  is  furthermore  shown  that  the  centrifugal 
acceleration  generated  by  rotations  of  the  principal  axes  of  the  Reynolds  stress  tensor  can 
make  the  second  derivative  singular  at  the  most  extreme  limits  of  realizable  turbulence.  This 
previously  overlooked  effect  appears  to  make  it  impossible  to  identically  satisfy  the  strong 
form  of  realizability  in  any  version  of  the  present  generation  of  second-order  closures.  On  the 
other  hand,  models  properly  formulated  to  satisfy  the  weak  form  of  realizability  -  wherein 
states  of  one  or  two  component  turbulence  are  not  accessible  in  finite  time  -  are  found  to 
be  realizable.  However,  unlike  the  simpler  and  more  commonly  used  second-order  closures, 
these  models  can  be  ill-behaved  near  the  extreme  limits  of  realizable  turbulence  due  to  the 
way  that  higher-degree  nonlinearities  are  often  unnecessarily  introduced  to  satisfy  realiz¬ 
ability.  Illustrative  computations  of  homogeneous  shear  flows  are  presented  to  demonstrate 
these  points  which  can  have  important  implications  for  turbulence  modeling. 
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1.  Introduction 


Second-order  closure  models  have  been  an  active  area  of  research  in  turbulence  model¬ 
ing  for  the  past  few  decades.  Since  second-order  closures  are  based  on  the  Reynolds  stress 
transport  equation  -  which  accounts  for  both  history  and  nonlocal  effects  -  these  models, 
in  principle,  allow  for  the  description  of  more  turbulence  physics  than  lower  level  closures. 
Schumann  (1977)  was  the  first  to  systematically  address  the  issue  of  realizability  in  second- 
order  closure  modeling.  The  constraint  of  realizability  requires  that  a  Reynolds  stress  model 
yield  non-negative  component  energies  in  all  turbulent  flows,  with  the  Schwarz  inequality 
satisfied  for  each  off-diagonal  component  of  the  Reynolds  stress  tensor.  Schumann  showed 
that  realizability  is  satisfied  identically  by  a  model  if,  starting  from  any  realizable  initial  con¬ 
ditions,  it  predicts  a  Reynolds  stress  tensor  with  non-negative  eigenvalues  for  all  later  times. 
He  also  provided  a  variety  of  necessary  conditions  for  the  satisfaction  of  realizability  and 
briefly  outlined  a  proposed  method  for  making  unrealizable  second-order  closures  realizable. 
This  issue  was  of  interest  since  it  had  long  been  known  that  unrealizable  second-order  clo¬ 
sure  models  can  lead  to  numerical  instabilities.  The  ad  hoc  numerical  technique  of  clipping 
-  whereby  when  a  negative  component  energy  is  computed  it  is  arbitrarily  set  to  zero  -  was 
introduced  to  alleviate  just  such  a  problem  (see  Deardorff  1973). 

Lumley  (1978,  1983)  was  the  first  to  advocate  the  systematic  use  of  realizability  con¬ 
straints  in  the  formulation  and  calibration  of  second-order  closure  models.  He  clarified  the 
constraints  that  a  second-order  closure  should  satisfy  in  order  to  be  consistent  with  the 
strong  form  of  realizability:  When  a  principal  Reynolds  stress  component  vanishes,  its  time 
rate  must  also  vanish  and  its  second  derivative  must  be  positive.  However,  Lumley  made 
a  significant  departure  from  Schumann  (1977)  in  so  far  as  he  suggested  that  realizability 
could  serve  as  a  powerful  new  constraint  in  determining  the  allowable  mathematical  form 
of  Reynolds  stress  models.  Lumley  (1978)  also  introduced  the  constraint  of  joint  realizabil¬ 
ity  into  second-order  closure  modeling  whereby  the  Schwarz  inequality  is  imposed  on  scalar 
fluxes.  A  few  years  later,  Shih  and  Lumley  (1985)  developed  a  rapid  pressure-strain  model 
based  on  the  implementation  of  invariant  tensor  theory  and  realizability  constraints  alone. 
They  combined  this  new  rapid  model  (after  some  subsequent  modifications  were  made)  with 
the  slow  pressure-strain  and  isotropic  dissipation  rate  models  that  Lumley  (1978)  had  de¬ 
veloped  earlier.  The  resulting  second-order  closure  has  been  commonly  referred  to  as  the 
Shih-Lumley  model  in  the  literature.  Shih  and  Lumley  -  who  claimed  that  their  model  was 
the  first  generally  realizable  second-order  closure  -  subsequently  reported  a  few  applications 
to  turbulent  shear  flows  (see  Shih  and  Lumley  1992,  1993).  However,  they  did  not  report 
any  tests  of  their  model  demonstrating  whether  it  did  indeed  guarantee  realizable  solutions. 

Pope  (1985)  departed  from  the  approach  of  Lumley  in  proposing  what  has  come  to  be 
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known  as  the  wecJc  form  of  reoHzabUity.  Pope  only  required  that  when  a  principal  Reynolds 
stress  component  vanishes,  its  time  derivative  be  positive.  This  does  guarantee  realizability 
in  that  non-negative  energy  components  cannot  occur  when  this  constraint  is  satisfied.  It 
also  has  the  advantage  of  eliminating  the  need  to  enforce  the  positive  second  time  derivative 
constraint,  since  a  model  that  properly  satisfies  the  weak  form  of  realizability  can  never  access 
one  or  two  component  states  of  turbulence  in  finite  time.  Lumley  and  co-workers  criticized 
the  weak  form  of  realizability  on  the  grounds  that  states  of  two-component  turbulence  -  which 
can  occur  in  practical  turbulent  flows  near  a  solid  boundary  -  were  inaccessible.  Nevertheless, 
Haworth  and  Pope  (1986)  and  Pope  (1993)  derived  second-order  closures  for  homogeneous 
turbulence  from  a  Langevin  equation  that  satisfied  this  weak  form  of  realizability  and,  hence, 
guaranteed  positive  component  enerfpes.  While  Langevin  equations  had  been  used  to  prove 
realizability  of  the  DIA  and  EDQNM  two-point  closures  (see  Kraichnan  1961  and  Orszag 
1970, 1977),  Pope  was  the  first  to  bring  this  stochastic  analysis  into  the  realm  of  second-order 
closure  modeling. 

The  purpose  of  the  present  paper  is  to  clarify  the  issue  of  realizability  in  second-order 
closure  modeling  and  to  discuss  alternative  means  of  ensuring  realizable  solutions.  In  regard 
to  the  former  point,  there  are  a  variety  of  confusing  and  conflicting  claims  in  the  literature 
that  need  to  be  clarified  -  a  task  that  can  be  accomplished  within  the  context  of  homogeneous 
turbulence.  It  will  be  proven  mathematically  and  demonstrated  computationally  that  the 
Shih-Lumley  model  yields  unrealizable  solutions.  The  lack  of  realizability  arises  from  the 
failure  to  satisfy  the  necessary  positive  second  time  derivative  constraint  when  a  principal 
Reynolds  stress  vanishes  -  a  crucial  condition  that  must  be  satisfied  by  models  that  allow 
access  to  one  or  two  component  states  of  turbulence.  This  problem  only  becomes  appsurent 
when  the  non-analytic  terms  in  the  Shih-Lumley  model  are  made  single- valued  as  required  for 
physical  consistency  (a  turbulence  model  cannot  contain  functions  that  are  multiple- valued 
and  give  rise  to  non-unique  solutions).  It  will  furthermore  be  shown  that  the  centrifugal 
acceleration  generated  by  rotations  of  the  principal  axes  of  the  Reynolds  stress  tensor  can 
make  second  and  higher-order  time  derivatives  singular  at  the  most  extreme  limits  of  two 
component  turbulence.  This  is  a  heretofore  neglected  effect  that  appears  to  make  any  version 
of  the  present  generation  of  second-order  closures  incapable  of  identically  satisfying  the  strong 
form  of  realizability.  Formulations  of  the  weak  form  of  realizability  -  which  do  give  rise  to 
realizable  models  -  are  examined  critically  and  a  more  physically  consistent  statement  of 
this  constraint  is  provided.  However,  models  that  satisfy  the  weak,  as  well  as  the  strong, 
form  of  realizability  can  be  ill-behaved  near  one  and  two  component  states  of  turbulence  due 
to  the  unnecessary  introduction  of  higher-degree  nonlinearities  to  satisfy  realizability.  In  a 
related  paper  (Durbin  and  Speziale  1993),  it  is  shown  how  linear  models  can  be  modified 


2 


via  a  stochastic  analysis  to  guarantee  realizability  in  such  extreme  cases  with  better  behaved 
predictions.  These  issues  will  be  discussed  in  detail  in  the  sections  to  follow  and  illustrative 
calculations  of  homogeneous  shear  flow  will  be  provided. 


2.  Realizability  Constraints  for  Second-Order  Closures 

We  will  consider  incompressible  turbulent  flows  governed  by  the  Navier-Stokes  and  con¬ 
tinuity  equations 

dvi  dVi  dP  _2 


where  v,-  is  the  velocity  vector,  P  is  the  kinematic  pressure,  and  v  is  the  kinematic  viscosity 
of  the  fluid.  As  in  all  studies  of  Reynolds  stress  modeling,  the  velocity  and  pressure  are 
decomposed  into  mean  and  fluctuating  parts  as  follows: 

Vi=Vi+Ui^  P  =  P  -f  p  (3) 

where  an  overbar  represents  an  ensemble  mean.  The  Reynolds-averaged  Navier-Stokes  and 
continuity  equations  take  the  form 

where  Tij  =  is  the  Reynolds  stress  tensor.  Prom  its  definition,  it  is  clear  that  r,j  has  non¬ 
negative  eigenvalues  -  a  property  that  realizability  constraints  seek  to  preserve  in  Reynolds 
stress  models. 

In  order  to  achieve  closure,  (4)-(5)  must  be  supplemented  with  a  Reynolds  stress  model 
that  ties  r,-j  to  the  global  history  of  the  mean  velocity  field  in  a  physically  reasonable  fashion. 
We  will  analyze  second-order  closure  models  that  are  based  on  the  Reynolds  stress  transport 
equation  which,  for  homogeneous  turbulence,  simplifies  to  the  form  (cf.  Hinze  1975) 


where 


/ dui  duj\  dui  duj 
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ate,  respectively,  the  pressure-strain  correlation  and  the  dissipation  rate  tensor.  The  dissi¬ 
pation  rate  tensor  can  be  split  into  isotropic  and  deviatoric  parts,  respectively,  as  follows: 


2  ^ 

£ij  —  3  D^ij 

(8) 

where 

_  dUi  dUi 

(9) 

^~‘'dXjdXj 

is  the  scalar  dissipation  rate.  Writing 

T,  2  ^ 

(10) 

with  n,j  =  —  D^ijt  it  follows  from  (6)  that  closure  is  achieved  once  models  for  Ilij  and  e 

are  provided.  The  current  generation  of  second-order  closures  are  based  on  models  for  11, j 
and  e  that  reduce  to  the  general  form  (cf.  Reynolds  1987  and  Speziale  1991) 


n.,  =  eAiiib)  -H  (11) 

e  =  (12) 

for  homogeneous  turbulent  flows.  In  (11)-(12), 


(13) 


are,  respectively,  the  turbulent  kinetic  energy  and  anisotropy  tensor;  the  coefficients  Cei  and 
Cc2  typically  taken  to  be  either  constants  or  functions  of  the  turbulence  Reynolds  number 
Ret  =  jve  and,  possibly,  the  dimensionless  invariants  of  dvijdzj. 

As  first  pointed  out  by  Schumann  (1977)  and  Lumley  (1978),  it  is  euier  to  examine 
the  question  of  realizability  in  a  coordinate  system  aligned  with  the  principal  axes  of  the 
Reynolds  stress  tensor.  However,  considerable  care  must  be  taken  in  such  an  analysis  since 
the  principal  axes  of  the  Reynolds  stress  tensor  can  rotate  in  a  time-dependent  manner  for 
temporally  evolving  homogeneous  turbulent  flows.  In  coordinate  free  notation,  (6)  can  be 
written  as  a  dynamical  system: 

f  =  /  (14) 


where 

f  =  p  +  n-\ti  (15) 

given  that  I  is  the  unit  tensor  and  P  is  the  production  tensor  whose  components  are  provided 
by  the  first  two  terms  on  the  right-hand-side  of  (6).  In  a  fixed  coordinate  system  with  base 


vectors  =  (ci,C2»C3),  the  components  of  (14)  are  given  by  (6).  However,  relative  to  the 
principal  ajces  -  which  have  the  base  vectors  =  (Ai,  A2,A3)  -  the  component  form  of  (14) 
is  more  complex  since  the  basis  is  rotating  in  time.  For  any  symmetric  second  rank  tensor 
T,  it  can  be  shown  that  (see  Appendix  A) 

{T)afi  =  Ta0  +  ea-,S^-,Ts0  +  eflyfil^Tsa  (16) 

where  ea0^  is  the  permutation  tensor  and  Hq  =  ^^(t)  is  the  angular  velocity  at  which  the 
principal  axes  are  rotating.  Relative  to  the  principal  axes,  (14)  takes  the  form 

~  fafi  CoTrS 

However,  the  Reynolds  stress  tensor  is  diagonal  relative  to  the  principal  axes: 

Ta0  (^®) 

where  r(aa)  (for  a  =  1,2,3)  are  the  principal  Reynolds  stresses  and  the  Einstein  summation 
convention  is  suspended  for  indices  that  lie  within  parentheses.  From  (17)  it  is  now  clear 
that 

^(aa)^a0  ~  fafi 

and,  thus,  we  have 

^(oa)  —  f(aa)  (20) 

by  setting  a  =  P  since  Ca-ya  =  0  from  the  definition  of  the  permutation  tensor.  It  has  thus 
been  shown  that,  consistent  with  the  earlier  proof  presented  by  Lumley  (1978),  rotations  of 
the  principal  axes  of  the  Reynolds  stress  tensor  have  no  effect  on  the  formulation  of  the  first 
derivative  constraint.  As  we  will  demonstrate  later,  such  rotations  have  an  important  effect 
on  the  formulation  of  the  second  derivative  constraint. 

Schumann  (1977)  showed  that  when  r(oa)  =  0,  it  follows  that  the  correlations 

P(aa)  =  0,  $(aa)  =  0,  £(00)  =  0  (21) 

as  a  rigorous  consequence  of  their  definitions  and  the  Schwarz  inequality.  Since 

f{aoi)  ~  P{ooi)  d”  (22) 


it  follows  that  f^aa)  vanishes  when  rjoa)  =  0.  This  leads  us  to  the  long  established  realizability 
constraint:  When 


^(aa)  —  0 


(23) 


it  follows  that 


T^aa)  =  0 


(24) 
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(see  Schumann  1977  and  Lumley  1978).  This  constraint  is  typically  satisfied  by  guaranteeing 
that  when  r(aa)  vanishes: 


,(5)  _ 

^(aa)  D^(qq)  — 


*(*)  _  0 
*(aa)  “ 


(25) 

(26) 


where  and  are,  respectively,  the  slow  and  rapid  parts  of  the  pressure-strain  correla¬ 
tion  (here,  #  =  -f  cf.  Lumley  1978).  While  the  realizability  constraints  (25)-(26) 
are  a  rigorous  consequence  of  the  Navier-Stokes  equations,  in  and  of  themselves,  they  are 
neither  necessary  nor  sufficient  to  guarantee  realizability.  As  shown  by  Pope  (1985),  realiz¬ 
ability  can  be  satisfied  if  (24)  is  replaced  with  the  weaker  constraint  f(aa)  '>  0.  Furthermore, 
even  if  the  constraint  f(aa)  =  0  is  satisfied  when  T(aa)  vanishes,  an  added  second  derivative 
condition  must  be  met  for  full  realizability;  namely,  we  must  have 


0  <  <  oo. 


(27) 


If  =  0,  then  the  first  non-vanishing  derivative  of  f(aa)  must  be  positive  and  bounded. 
It  is  crucial  that  the  non-negative  second  derivative  constraint  (27)  be  properly  satisfied  to 
guarantee  the  realizability  of  second-order  closures  that  allow  one  or  two  component  states  of 
turbulence  to  be  accessible  in  finite  time  (i.e.,  (27)  is  a  critical  constraint  for  the  satisfaction 
of  the  strong  form  of  realizability).  This  constraint  has  not  been  properly  analyzed  in  previous 
applications  of  realizability  to  second-order  closure  modeling,  as  we  will  show  in  the  next 
section. 


3.  Analysis  of  Existing  Second-Order  Closures 

In  this  section,  the  strong  form  of  realizability  will  be  critically  assessed  based  on  an 
analysis  of  the  model  proposed  by  Shih  and  Lumley.  A  new  statement  of  the  weak  form  of 
realizability  will  also  be  provided  along  with  am  example  of  a  model  (the  Fu,  Launder  and 
Tselepedakis  model)  that  satisfies  this  constraint  and  is  realizable. 


3.1.  The  Strong  Form  of  Realizability 


In  the  Shih-Lumley  model  (see  Appendix  B),  an  attempt  is  made  to  satisfy  the  second 
derivative  constraint  (27)  by  the  construction  of  a  rapid  pressure-strain  model  that  behaves 
asymptotically  as 


.  1/2 
T(aa)  OC  T^aa) 


(28) 


near  two-component  states  of  turbulence  where  r^aa)  ^  0.  Here,  it  is  understood  that  the 
proportionality  factor  must  be  negative  in  order  to  aUow  access  to  two-component  states. 
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Tensorially  invariant  forms  that  behave  like  (28)  near  T(oq)  w  0  (where  a  =  1,2  or  3)  were 
obtained  by  Shih  and  Lumley  by  the  introduction  of  the  invariant  function 

F  =  1  +  911  +  271 1 1  (29) 

where  JI  and  III  are,  respectively,  the  second  and  third  invariants  of  the  anisotropy  tensor 
bij  (see  Appendix  B  and  Lumley  1978).  It  is  a  simple  matter  to  show  that 

^  glnn.mn.33i 

in  terms  of  the  principal  Reynolds  stresses.  For  any  realizable  turbulence,  we  must  have 
0  <  F  <  1  (see  Lumley  1978).  This  prompted  Shih  and  Lumley  (1985)  to  suggest  the 
tensorially  invariant  form 

ncc)  OC  (31) 

in  the  limit  as  T(Qa))  and  hence  F,  goes  to  zero.  The  full  tensorial  asymptotic  form  of  the 
Shih-Lumley  model  near  =  0  is  given  by 

(is  ‘ 

where  V  =  —Tijdvifdxj  is  the  turbulence  production  and  ai  and  are  model  constamts 
that  are  positive.  Since  the  coefficient  that  multiplies  F^!^  in  (32)  can  be  negative,  the 
Shih-Lumley  model  allows  access  to  one  and  two-component  states  of  turbulence. 

The  One-Component  Limit 

The  logic  used  by  Shih  and  Lumley  appears  on  the  surface  to  be  sound.  They  argued 
that  (28)  guarantees  that  when  r(oQ)  vanishes,  it  follows  that  T(aQ)  vanishes  with  T(aa)  >  0  - 
a  set  of  conditions  which  shoiild  ensure  that  T^aa)  never  becomes  negative  as  discussed  in  the 
previous  section.  While  this  line  of  reaisoning  appears  to  be  correct  at  first  glance,  a  deeper 
analysis  uncovers  a  fundamental  problem.  By  making  use  of  (30),  it  is  clear  that  for  r(aa) 
and  F  close  to  zero,  (32)  can  be  rewritten  in  the  form 

F  =  -CF^/*  (33) 

where  in  the  neighborhood  of  r(aa)  =  0  we  can  treat  the  coefficient  (7  as  a  constant.  On  the 
basis  of  (33)  it  follows  that  when  F  =  0:  F  =  0  and  F  =  which  is  positive.  However,  as 
we  will  soon  see,  it  would  be  a  mistake  to  conclude  that  (33)  guarantees  realizability  due  to  a 
variety  of  difficulties.  Even  though  (33)  can  yield  solutions  where  F  >  0  for  all  times,  there 
is  a  subtle  problem.  While  F  <  0  corresponds  to  unrealizable  turbulence,  the  converse  is  not 
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true  for  ^  >  0  (see  Figure  1).  A  turbulence  that  becomes  unrealizable  by  passing  through 
a  one-component  state  can  correspond  to  F  >  0.  For  example,  consider  a  first  order  Taylor 
expansion  for  b^j  in  the  neighborhood  of  a  one  component  turbulence  where  rn  =  r22  =  0: 


/-l-at 


0 

0 


0  0 
0 

0  I  -I-  (a  -I-  /S)t  J 


(34) 


Here  a  and  ^  are  constants  and  t  is  a  short  elapsed  time  after  a  one-component  state  is 
achieved.  When  either  a  ot  w  greater  than  zero,  the  turbulence  becomes  unrealizable 
since  —  ^  ^  &(»)  <  |  for  realizable  turbulence.  It  is  a  simple  matter  to  show  from  Eq.  (30) 
that  if  a  >  0  and  /3  =  0,  then  F  =  0;  furthermore,  when  a  >  0  and  >  0,  then  F  >  0. 
In  either  case,  the  turbulence  becomes  unrealizable  with  F  >  0.  This  leads  us  to  our  first 
pertinent  conclusion:  The  methodology  of  Shih  and  Lumley  is  fundamentally  incapable  of 
guaranteeing  realizability  for  turbulent  flows  that  are  near  a  one-component  state.  While 
this  is  not  the  major  deficiency  with  the  Shih-Lumley  approach,  it  is  still  of  consequence 
since  turbulence  that  is  near  a  one-component  state  is  asymptotically  approached  in  certain 
geophysical  flows  with  strong  stratification  (see  Zeman  and  Lumley  1976).  A  model  that 
claims  to  be  fuUy  realizable  must  accommodate  such  a  limit. 


The  Central  Problem  of  Non-Analyticity 

The  major  deficiency  with  the  Shih-Lumley  model  lies  in  its  use  of  the  non-analytic 
function  F^^^  as  a  means  to  satisfy  realizability.  If  the  dimensionless  time  t*  =  J  Cdt  is 
introduced  into  (33),  it  follows  that  its  exact  solution  can  be  written  in  the  form 

F^/2  =  (35) 

Since  F^/^  =  ±\/F  (where  yfF  denotes  the  positive  value  of  the  square  root,  i.e.  y/F  = 
|F^/^|),  it  follows  that  there  are  two  principal  brsmches  to  the  solution: 

=  (\/ft  +  (36) 

and 

F  =  Iff  (37) 

which  are  illustrated  graphically  in  Figure  2  for  an  initial  condition  of  Fo  =  1.  If,  starting 
from  any  Fq  >  0,  we  want  to  pass  through  a  two-component  state  of  turbulence  (where 
F  =  0)  and  then  turn  back  up,  we  must  pick  the  second  branch  of  this  solution  given  by 
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(37).  However,  there  is  a  problem:  for  0  <  <  2,  we  have  >  0  whereas  for  >  2  we 

have  <  0  as  a  result  of  (35).  Hence,  it  is  clear  that  the  solutions  of  (33)  are  remaining 
realizable  by  virtue  of  shifting  branches  between  y/F  and  —\/F.  It  is  unacceptable  for 
a  turbulence  model  to  have  multiple  valued  functions;  when  calculating  a  complex  turbulent 
flow,  how  do  we  know  when  to  take  F^^^  equal  to  \/F  or  — \/F?  There  is  no  physically  based 
selection  rule.  The  differential  equations  corresponding  to  the  two  branches  of  F^^^  are 

F  =  \/F  (38) 

and 

F  =  -y/F  (39) 

Eq.  (39)  has  no  solution  for  t'  >  2  (the  solutions  to  (38)  and  (39)  for  0  <  t*  <  2  are 
illustrated  with  the  solid  lines  in  Figure  2).  Shih  and  Lumley  tacitly  take  F*'^^  =  y/F  (if 
they  were  to  take  F^^^  =  —y/F,  then  two  component  states  would  not  be  accessible  when 
C  >  0]  furthermore,  their  model  then  predicts  equilibrium  values  of  the  normal  Reynolds 
stress  anisotropies  that  are  of  the  wrong  sign).  Hence,  the  Shih-Lumley  model  actually 
behaves  like  the  generic  differential  equation  (39)  in  the  vicinity  of  F  =  0  (for  C  >  0)  which 
has  no  solution  when  t*  >  2.  A  differential  equation  like  (39)  that  yields  an  undefined  F  is 
technically  unrealizable  since  we  must  have  0  <  F  <  1  in  a  realizable  turbulence. 

The  fact  that  there  is  a  realizability  problem  with  (39)  can  more  e^lsily  be  seen  in  an 
alternative  way.  If  (39)  had  realizable  solutions  for  all  times  where  F  >  0,  then  we  should 
be  able  to  replace  y/F  in  (39)  with  the  ^|F|.  The  stable  solution  to  the  differential  equation 

F  =  (40) 

is  illustrated  in  Figure  3  and  it  is  decidedly  unrealizable:  F  — >  — oo  as  t*  oo  (although 
(40)  formally  has  the  fixed  point  F  =  0,  it  is  unstable).  The  results  displayed  in  Figure  3 
were  obtained  by  a  Runge-Kutta  numerical  integration.  It  can  be  shown  that  the  analytical 
form  of  the  stable  solution  to  (40)  is  given  by 

F  =  (I  -  2//(f  -  2)]  (41) 

where  //(•)  is  the  Heaviside  function.  From  (41)  it  follows  that  F  is  discontinuous  at  the 
critical  time  t*  =  2  when  F  vanishes;  the  one-sided  second  derivatives  are  F~  =  ^  and 
F"*"  =  —\  when  F  =  0.  For  the  nonlinear  differential  equation  (39),  F~  =  |  whereas  F+ 
does  not  exist  ai  t*  —  2  where  F  =  0.  In  either  case,  F  is  undefined  when  F  =  0  and, 
hence,  the  crucial  second  derivative  constraint  (27)  is  not  satisfied.  This  leads  us  to  the 
following  important  conclusion:  When  the  non-analytic  terms  in  the  Shih-Lumley  model  are 
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made  single-valued,  as  required  on  physical  grounds,  it  fails  to  satisfy  the  crucial  positive 
second  derivative  condition  -  a  fatal  flaw  that  can  lead  to  realizability  violations.  In  the 
next  section,  illustrative  computations  of  homogeneous  shear  flow  will  be  presented  which 
definitively  demonstrate  that  the  Shih-Lumley  model  is  not  realizable. 


Rotations  of  the  Principal  Axes 

The  use  of  the  non-analytic  term  is  not  a  viable  means  for  ensuring  realizability  as 
shown  above.  The  question  remains  as  to  whether  there  exist  alternative  mathematical  forms 
that  can  identically  satisfy  the  strong  form  of  realizability.  It  will  now  be  shown  that  the 
centrifugal  acceleration  generated  by  rotations  of  the  principal  axes  of  the  Reynolds  stress 
tensor  can  make  T(aa)  singular  at  two  of  the  most  extreme  limits  of  realizable  turbulence  in 
any  version  of  the  present  generation  of  second-order  closures.  If  f(aa)  is  singular,  realizability 
cannot  be  guaranteed,  ('or  example,  consider  the  differential  equation 

F  =  (42) 

Eq.  (42)  has  the  exact  solution 

(«) 

which  becomes  negative  and  violates  realizability.  It  is  a  simple  matter  with  Eq.  (42)  to 
show  that  when  F  =  0:  F  =  0  and  F  =  -i-oo.  Hence,  a  singular  second  derivative,  even  if  it 
is  positive,  can  be  fatal  with  respect  to  realizability. 

By  differentiating  (14),  we  obtain  the  equation 

f  =  /.  (44) 


From  Appendix  A,  it  follows  that  the  component  form  of  (44)  relative  to  the  principal  axes 
is  given  by: 

Taff  —  2ep.ygil.yTga 


(45) 


-i-2(n.,n^T.yo-  n*7Va-)^a^  900 

•  • 
where  gap  =  (/)aj9  (we  write  this  equation  in  terms  of  gap  since,  unlike  fap,  it  is  free  of  any 

explicit  dependence  on  fJ).  Since,  relative  to  the  principal  axes,  we  have  Tap  =  T{aa)^Q0  and 


10 


i'aff  =  T(ao,)Sa0i  it  foUows  that 

'5^(aa)  ~  +  4n*T(aa)  +  2(trT)fiJ„) 


+2 


(trr) 


(46) 


■f  9(aa) 


where  and  trr  =  rjnj  +  T(22)  +  '’'(aa)-  *ii  existing  second-order 

closure  models,  fij  is  taken  to  be  of  the  form: 


/«  =  /j(’-,S,W,£) 


(47) 


where,  relative  to  the  fixed  coordinate  system  Zi,  the  components  of  S  and  are  given  by 

-q  =  i  4.  -  i  f  ^ 

2  \dxj  dxi)  ’  2  Uzj  azj  ■ 

Here,  we  are  considering  homogeneous  turbulent  flows  where  Sij  and  uJij  are  constant  tensors. 
However,  while  Sij  and  a7,j  are  zero,  the  same  is  not  true  for  Sa0  aaid  HJaff  which  are  non-zero 
and  depend  on  f)  due  to  the  rotation  of  the  principal  axes.  Hence,  depends  explicitly 
on  n  with  a  functional  form  determined  by  (47)  and  the  homogeneous  turbulence  under 
consideration.  In  order  to  isolate  the  terms  arising  from  the  rotation  of  the  principal  axes, 
we  wrote  (45)  in  terms  of  ga^  since: 


after  (44)  and  (47)  are  implemented.  It  is  clear  from  (49)  that  g^aa)  does  not  depend  explicitly 
on  n  (this  stands  in  contrast  to  f{aa)  w^hich  depends  on  Sap,  u’a/J  ^d,  hence,  on  as  a 
result  of  (16)).  The  off-diagonal  components  of  (19)  give 


■^(22)  ~  ^(33)  ■’^(33)  ~  ’■(11) 


11(3)  = 


-  /(»2) 
’•(11)  -  r(22) 


(50) 


Hence,  fl  can  become  singular  when  a  principal  component  of  the  Reynolds  stress  tensor, 
say  ’■(11),  vanishes;  this  will  happen  when  r(22)  or  r(33)  also  vanishes  (the  one-component 
limit)  and  when  r(22)  =  ’’(33)  (the  axisymmetric  two-component  limit).  At  these  two  critical 
points  -  which  constitute  the  end  points  of  the  two  component  line  of  the  Lumley  anisotropy 
invariant  map  -  the  second  derivative  f(oQ)  can  become  singular  according  to  (46)  and  (50). 
This  makes  it  impossible  to  satisfy  th'2  crucial  constraint  (27). 
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We  thus  conclude  that  it  appears  to  be  impossible  to  identically  satisfy  the  strong  form 
of  realizability  in  the  current  generation  of  second-order  closure  models.  The  only  way  to 
unequivocally  guarantee  realizability  in  the  current  version  of  second-order  closures  is  to 
avoid  access  to  one  or  two  component  states  of  turbulence  in  the  homogeneous  limit. 

3.2.  The  Weak  Form  of  Realizability 

We  will  now  consider  the  weak  form  of  realizability  where  one  or  two  component  states 
are  inaccessible  in  finite  time.  As  first  presented  by  Pope  (1985),  weak  realizability  is  satisfied 
if,  when  a  principal  Reynolds  stress  r^aa)  vanishes, 

T(c,a)  >  0.  (51) 

This  condition  ensures  that  an  initially  three-component  turbulence  never  achieves  a  two- 
component  state.  Realizability  is  guaranteed  without  the  need  to  impose  a  second  derivative 
constraint.  Haworth  and  Pope  (1986)  derived  a  second-order  closure  based  on  a  Langevin 
equation.  Weak  realizability  is  satisfied  therein  as  follows:  When  a  principal  Reynolds  stress 
Tf^aa)  vanishes, 

=  0  (52) 

♦if!)  -  >  O-  (53) 

In  the  Haworth  and  Pope  (1986)  model,  oCij  =  0,  so  it  follows  that  (52)-(53)  ensures 
the  satisfaction  of  the  weak  realizability  constraint  (51).  In  that  model  -  as  in  virtually 
every  version  of  the  current  generation  of  second-order  closures  -  the  slow  pressure  strain 
correlation  is  modeled  as 

=  -C.ebij  -H  C2e{bi,bkj  -  ^bMij)  (54) 

where  the  coefficients  Ci  and  C2  can  be  functions  of  the  invariants  II  and  III  as  well  as  the 
turbulence  Reynolds  number.  Since  b^aa)  =  —5  when  r(o„)  =  0,  the  constraint  (53)  can  be 
satisfied  when  a  principal  Reynolds  stress  component  vanishes  provided  that  (see  Sarkar  and 
Speziale  1990) 

Cl  >2,  Cj  <  3(C',  -  2).  (55) 

Hence,  any  slow  pressure-strain  model  of  form  (54)-(55),  when  combined  with  a  rapid 
pressure-strain  model  that  vanishes  in  the  two  component  limit,  will  be  realizable. 

Lumley  and  co-workers  have  criticized  weak  realizability  on  the  grounds  that  one  and 
two  component  states  are  made  inaccessible;  they  often  cite  the  two-component  turbulence 
that  occurs  near  a  solid  boundary  as  a  counter-example.  We  do  not  consider  this  criticism  to 
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be  compelling.  Weak  realizability  only  precludes  an  initially  three-component  homogeneous 
turbulence  £rom  achieving  a  one  or  two  component  state  in  finite  time;  it  renders  no  constraint 
on  inhomogeneous  turbulent  flows  such  as  those  which  occur  near  a  solid  boundary.  The 
more  serious  criticism  of  weak  realizability  is  that  (51)  is  inconsistent  with  the  Navier-Stokes 
equations.  As  proven  by  Schumann  (1977),  when  a  principal  Reynolds  stress  component 
r(Q,Q)  vanishes,  we  must  have  T(aa)  =  0.  This  can  be  easily  seen  from  the  associated  Schwarz 
inequality: 

T(aa)  =  (56) 

From  (56)  it  is  clear  that  T{aa)  vanishes  when  r^aa)  vanishes.  Hence,  we  propose  an  alternative 
form  of  weak  realizability:  When  a  principal  Reynolds  stress  ^aa)  vanishes,  we  require  that 
■^(aa)  =  0>  however,  when  is  in  a  neighborhood  of  zero,  we  then  enforce  the  constraint 
that  ‘r(aa)  >  0.  The  latter  condition  guarantees  that  an  initially  three-component  turbulence 
never  achieves  a  one  or  two  component  state;  however,  it  allows  model  predictions  to  come 
arbitrarily  close  to  the  boundaries  of  realizable  turbulence  unlike  the  formulation  of  Pope 
(1985).  This  alternative  form  of  weak  realizability  might  be  satisfied  as  follows:  when  r^aa) 
is  arbitrarily  close  to  zero,  enforce  the  constraints 

«lfl)  =  f  (57) 

=  0(F*)  (68) 

with  b  >  a  and  C  >  0.  In  traditional  slow  pressure-strain  models,  (57)  corresponds  to  the 
conditions 

Cl  =  2  +  0(f’“)  >  2  (59) 

C2  <  3(C'i  -  2).  (60) 

Here,  the  exponents  a  and  6  are  completely  arbitrary  (for  simplicity,  we  have  absorbed  peij 
into  Interestingly  enough,  the  model  of  Fu,  Launder,  and  Tselepidakis  (1987)  (see 

Appendix  B)  satisfies  this  alternative  form  of  weak  realizability  with  the  exponent  a  =  | 
and  6=1.  In  the  next  section,  we  will  demonstrate  computationally  that  this  model  is 
realizable,  i.e.,  it  never  3rields  negative  component  energies. 

4.  Illustrative  Computations 

The  theoretical  points  discussed  in  the  previous  section  will  now  be  demonstrated  com¬ 
putationally  for  homogeneous  shear  flow.  The  test  case  of  homogeneous  shear  flow  is  selected 
since  it  constitutes  one  of  the  simplest  and  most  important  benchmark  turbulent  flows  that 
has  been  documented  extensively  by  physical  and  numerical  experiments  (  Tavoularis  and 
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Corrsin  1981,  Tkvoularis  &nd  Karnik  1989,  and  Rogers,  Moin  and  Reynolds  1986).  In  these 
experiments,  an  initially  isotropic  turbulence  evolves  after  it  is  subjected  to  a  uniform  shear 
rate  5.  However,  in  order  to  provide  a  stringent  test  of  realizability,  we  will  consider  initial 
conditions  that  are  anisotropic  and  close  to  a  one  or  two  component  state  of  turbulence 
where  F  m  0  (the  existing  physical  and  numerical  experiments  as  mentioned  above  have 
concentrated  on  isotropic  initial  conditions  where  F  =  1).  Furthermore,  we  will  consider 
large  initial  values  of  the  shear  parameter  SK/e  which  can  make  a  model  more  prone  to 
experience  realizability  violations. 

The  computations  of  homogeneous  shear  flow  to  be  presented  consist  of  numerical  solu¬ 
tions  of  the  nonlinear  ordinary  differential  equations  (6)  and  (12)  subject  to  the  constant 
mean  velodty  gradient  tensor 

P-  =  SSi^6i2  (61) 

CfXj 

and  the  initial  conditions 

K  =  (6.i)o.  ^  ^  (62) 

€  €o 

at  time  t  =  0.  Three  models  will  be  considered  for  n,’j,  Cei  and  0^2  in  Eqs.  (6)  and  (12):  the 
Shih-Lumley  (SL)  model;  the  Fu,  Launder  and  Tselepidakis  (FLT)  model;  aind  the  IP  model 
which  is  a  simplified  form  of  the  Launder,  Reece  and  Rodi  (1975)  model  that  has  been 
used  in  a  variety  of  applications  starting  with  Gibson  and  Launder  (1978).  The  detailed 
form  of  these  models  is  provided  in  Appendix  B.  Time  accurate  solutions  were  computed 
using  a  fourth-order  accurate  Runge-Kutta  numerical  integration  scheme.  As  in  all  existing 
second-order  closures,  the  models  predict  that  bij  and  SK/e  eventually  achieve  equilibrium 
values  that  are  independent  of  the  initial  conditions  (these  constitute  the  fixed  points  of  the 
nonlinear  differential  equations  (6)  and  (12)).  Each  of  the  models  yield  realizable  fixed  points 
for  homoge"^us  shear  flow.  A  comparison  of  the  equilibrium  predictions  of  each  of  these 
three  models  with  physical  and  numerical  experiments  can  be  found  in  Speziale,  Gatski  and 
Mac  GioUa  Mhuiris  (1990)  and  Speziale,  Gatski  and  Sarkar  (1992). 

First,  we  will  present  computed  results  for  the  initial  conditions 

(jK 

ibn)o  =  (622)0  =  -0.32,  - ^  =  50  (63) 

Co 

given  that  (5,j)o  =  0  for  i  ^  y  (here,  and  in  the  other  computations  to  follow,  the  corre¬ 
sponding  initial  condition  for  633  can  be  obtained  from  the  traceless  constraint  6,-,-  =  0).  This 
corresponds  to  a  strong  uniform  shear  applied  to  an  initially  anisotropic  turbulence  that  is 
close  to  a  one-component  state  (the  reader  should  remember  that  if  611  —  622  =  —1/3,  then 
Til  =  T22  =  0).  Figure  4(a)  shows  the  computed  tr^i  jectories  of  the  Shih-Lumley  model  in  the 
phase  space  (—//,///)  corresponding  to  these  initial  conditions  for  homogeneous  shear  flow. 
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The  Lumley  ftnisotropy  invariant  map  is  superimposed  on  this  figure;  realizable  turbulence 
lies  within  this  curvilinear  triangle  where  the  origin  corresponds  to  isotropic  turbulence, 
point  A  corresponds  to  two-component  axisymmetric  turbulence,  and  point  B  corresponds 
to  one-component  turbulence  (see  Lumley  1978).  It  is  clear  from  Figure  4(a)  that  the  Shih- 
Lumley  model  becomes  unrealizable  since  its  solution  exits  the  triangle.  The  computations 
were  conducted  by  setting  in  the  rapid  model;  this  is  the  only  way  that  their 

model  is  both  single-valued  and  computable  to  an  equilibrium  state.  As  alluded  to  earlier, 
if  we  set  =  y/F^  their  model  has  no  solution  after  F  passes  through  zero  (when  solved 
numerically  vdth  any  finite  time  step,  a  negative  value  of  F  is  eventually  computed  before 
the  critical  time  where  F  vanishes  is  arrived  at  -  an  occurrence  that  causes  the  program  to 
terminate  in  an  error).  The  computed  results  were  validated  numerically  by  an  exhaustive 
grid  refinement  study. 

In  contrast  to  these  results,  the  IP  model  of  Launder  and  co-workers  is  fully  realizable  for 
this  set  of  initial  conditions  as  shown  in  Figure  4(b).  The  time  evolution  of  F  predicted  by  the 
Shih-Lumley  and  IP  models  is  compared  in  Figure  5(a)  where,  henceforth,  the  dimensionless 
time  t*  =  St.  Two  noteworthy  conclusions  can  be  drawn  from  these  results:  (a)  The  Shih- 
Lumley  model  is  unrealizable  (F  <  0)  for  the  large  elapsed  time  of  0  <  <  15,  and  (b) 

the  Shih-Lumley  model  undergoes  large  amplitude  oscillations  until  St  120  -  a  result  that 
will  be  shown  later  to  be  unphysical.  Similar  conclusions  can  be  drawn  from  Figure  5(b)  for 
the  time  evolution  of  the  shear  anisotropy  612.  The  Shih-Lumley  model  oscillates  between 
positive  and  negative  values  which  appears  to  be  highly  unphysical  (positive  values  of  612 
correspond  to  negative  turbulence  production). 

One  might  be  tempted  to  argue  that  the  initial  conditions  (63)  are  extreme  and,  therefore, 
constitute  an  overly  stringent  test  of  the  Shih-Lumley  model.  However,  it  should  be  said  at 
the  outset  that  a  model  which  claims  to  satisfy  the  strong  form  of  realizability  must  withstand 
such  a  test  -  particularly,  since  the  ol'der  IP  model,  which  makes  no  claim  of  being  generally 
realizable,  is  able  to  do  so.  Furthermore,  even  if  we  relax  the  initial  conditions  (63)  somewhat 
to  the  alternative  values: 

SK 

(tii)o  =  —0.3,  (622)0  =  —0.3,  - =  10  (64) 

Co 

the  Shih-Lumley  model  still  predicts  unrealizable  results  as  shown  in  Figure  6(a).  The  model 
remains  ill-behaved,  yielding  large  amplitude  oscillations  for  the  time  evolutions  of  F  and 
612,  as  shown  in  Figures  6(b)  -  6(c).  We  also  conducted  a  variety  of  calculations  for  initial 
conditions  that  were  near  the  center  of  the  t  cotroonent  line  of  the  Lumley  anisotropy 
invariant  map.  In  Figure  7(a),  the  trajectories  in  the  phase  space  {—II,  III)  predicted  by 
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the  Shih-Lumley  model  axe  shown  for  the  initial  conditions 

SK 

(611)0  =  -0.32,  (622)0  =  -0.12,  — ^  =  50  (65) 

Co 

in  homogeneous  shear  flow  (again,  (6,-j)o  =  0  for  t  ^  j).  The  Shih-Lumley  model  becomes 
unrealizable  whereas  for  these  same  initial  conditions,  the  IP  model  again  yields  realizable 
results  as  shown  in  Figure  7(b).  Other  computations  indicated  that  the  Shih-Lumley  model 
can  be  driven  unrealizable  for  virtually  any  initial  conditions  that  are  near  a  two-component 
state  of  turbulence  -  results  that  are  consistent  with  the  analysis  presented  in  the  previous 
section. 

As  mentioned  earlier,  the  IP  model  is  not  a  fully  realizable  second-order  closure.  However, 
the  only  way  that  the  IP  model  can  be  made  to  yield  unrealizable  results  in  homogeneous 
shear  flow  is  for  the  initial  production  to  be  negative  and  of  a  sufficient  magnitude.  This  is 
easy  to  see  once  we  recognize  that  when  a  principal  Reynolds  stress  r^aa)  vanishes  the  IP 
model  yields  (see  Appendix  B) 

n(„.)  -  =  it7.£  +  (66) 

which  become  negative  only  if  P  <  (l  —  e  (this  then  yields  f(o,a)  <  0  which  leads 
to  realizability  violations).  In  Figure  8(a)  the  trajectories  in  the  phase  space  (-//,///) 
obtained  from  the  IP  model  are  shown  for  the  initial  conditions: 

SK 

{bn)o  =  -0.24,  (622)0  =  0.17,  (612)0  =  0.2,  - -  =  50.  (67) 

eo 

Since  (612)0  is  positive,  the  initial  production  is  negative  and  the  IP  model  becomes  unre¬ 
alizable.  On  the  other  hand,  the  FLT  model  -  which  is  a  nonlinear  extension  of  the  IP 
model  that  was  formulated  to  satisfy  realizability  constraints  -  yields  realizable  results  for 
these  initial  conditions  as  shown  in  Figure  8(b).  Although  it  is  realizable,  the  FLT  model 
has  problems  with  large  amplitude  oscillations  like  those  of  the  Shih-Lumley  model.  This  is 
illustrated  in  Figures  9(a)-9(b)  where  the  predictions  of  the  FLT  model  for  F  and  612  are 
compared  with  results  obtained  from  the  IP  model.  The  IP  model  is  only  unrealizable  for 
a  relatively  short  transient  (i.e.,  for  0  <  <  2).  Although  the  FLT  model  is  realizable,  its 

solution  is  contaminated  by  large  amplitude  oscillations  until  St  »  15. 

In  order  to  further  illustrate  the  problem  that  nonlinear  models  such  as  the  Shih-Lumley 
and  FLT  models  have  with  large  amplitude  oscillations,  we  return  to  the  near  one-component 
initial  conditions  (63).  For  this  case,  the  FLT  model  is  realizable  as  demonstrated  by  the 
phase  space  trajectories  shown  in  Figure  10.  However,  some  interesting  conclusions  can  be 
drawn  from  a  comparison  of  the  ratio  of  production  to  dissipation  {Vie)  predicted  by  the 
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IP,  FLT  and  Shih-Lumley  models  (see  Figures  ll(a)-(c)).  The  IP  model  predicts  a  "Pje  that 
peaks  at  5t  10  and  then  monotonically  approaches  its  equilibrium  value  oiVJe  »  2  (see 
Figure  11(a)).  While  we  have  no  confidence  in  the  quantitative  accuracy  of  these  results, 
the  crucial  point  is  that  the  IP  model  predicts  a  strictly  non-negative  turbulence  production. 
On  the  other  hand,  the  Shih-Lumley  and  FLT  models  predict  that  V/e  undergoes  large 
amplitude  oscillations  between  positive  and  negative  values  (see  Figures  11(b)  -  11(c)).  We 
believe  that  these  results  are  highly  unphysical;  there  is  no  apparent  physical  mechanism 
by  which  a  homogeneous  shear  flow  can  produce  a  negative  turbulence  production,  given 
that  it  is  initially  non-negative.  Further  evidence  concerning  the  unphysical  nature  of  large 
amplitude  oscillations  in  homogeneous  shear  flow  can  be  obtained  &om  Rapid  Distortion 
Theory  (RDT)  which  we  would  expect  to  be  an  excellent  approximation  to  the  Navier- 
Stokes  equations  for  homogeneous  shear  flow  when  SKofco  =  50  (see  Lee,  Kim  and  Moin 
1990).  In  Figure  12,  the  prediction  of  the  FLT  model  for  the  time  evolution  of  r22/(r22)o  is 
compared  vdth  the  RDT  solution  for  the  initial  conditions: 

SK 

(ill  )o  = -0.0833,  (622)0  = -0.0833,  — ^  =  50  (68) 

£0 

where  (6,^)0  =  0  for  t  ^  j.  These  initial  conditions  correspond  to  an  axisymmetric  homoge¬ 
neous  turbulence  for  which  we  were  able  to  compute  an  energy  spectrum  tensor  -  an  input 
that  is  needed  as  an  initial  condition  for  RDT.  It  is  clear  from  the  results  in  Figure  12  that 
the  FLT  model  is  predicting  spurious  osdllations  (the  same  is  true  of  the  Shih-Lumley  model 
which  is  not  shown).  These  spurious  oscillations  appear  to  arise  from  the  introduction  of 
higher-degree  nonlinearities  in  6,^  within  the  rapid  pressure-strain  models  (see  Appendix  B). 
It  must  be  remembered  that  the  rapid  pressure-strain  correlation,  from  its  definition,  is  a  lin¬ 
ear  functional  of  the  energy  spectrum  tensor.  The  IP  model  -  which  unlike  the  Shih-Lumley 
and  FLT  models  is  consistent  with  this  linear  property  -  does  not  experience  unphysical 
oscillations.  Evidence  is  beginning  to  accumulate  concerning  the  counter-productive  effect 
that  higher  degree  nonlinearities  have  on  the  performance  of  pressure-strain  models. 

5.  Conclusions 

The  results  of  a  detailed  study  on  the  realizability  of  second-order  closure  models  has 
been  presented  within  the  framework  of  homogeneous  turbulence.  Several  interesting  and 
surprising  conclusions  were  arrived  at  which  can  be  summarized  as  follows: 

(1)  The  Shih-Lumley  model  is  not  a  realizable  model  and,  in  many  instances,  suffers  more 
realizability  violations  than  the  older  and  more  widely  used  IP  model  of  Launder  and 
co-workers.  This  lack  of  realizability  becomes  apparent  when  the  non-analytic  terms  in 
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their  model  containing  are  made  single-valued  in  a  physically  consistent  fashion 
-  a  simplification  which  invariably  leads  to  a  violation  of  the  crucial  positive  second 
derivative  constraint  which  requires  that  0  <  r^aa)  <  oo  when  T^aa)  vanishes.  Since  the 
Shih-Lumley  model  was  calibrated  largely  based  on  what  now  appear  to  be  questionable 
realizability  considerations,  it  becomes  more  understandable  why  the  model  does  not 
perform  well  in  basic  benchmark  turbulent  flows  as  recently  shown  by  Speziale,  Gatsld 
and  Sarkar  (1992)  and  Abid  and  Speziale  (1993). 

(2)  The  centrifugal  acceleration  generated  by  rotations  of  the  principal  axes  of  the  Rejmolds 
stress  tensor  can  render  singular  at  the  limits  of  one  component  turbulence  and 
axisymmetric  two  component  turbulence  (the  endpoints  A  and  B  of  the  two-component 
line  of  the  Lumley  anisotropy  invariant  map  shown  in  Figure  4(a)).  This  makes  it 
impossible  to  guarantee  the  satisfaction  of  the  positive  second  derivative  constraint 
(27)  in  existing  second-order  closure  models  -  a  deficiency  that  can  lead  to  realizability 
violations.  It  thus  appears  to  be  impossible  to  identically  satisfy  the  strong  form  of 
realizability  in  the  current  generation  of  second-order  closures. 

(3)  The  only  way  to  unequivocally  guarantee  realizability  in  the  current  generation  of 
second  order  closures  is  via  the  weak  form  of  realizability  where  access  to  one  or 
two  component  states  of  turbulence  is  avoided.  The  weak  form  of  realizability  as 
first  proposed  by  Pope  (1985)  does  indeed  guarantee  non-negative  component  energies 
in  arbitrary  homogeneous  turbulent  flows.  However,  there  is  an  associated  Schwarz 
inequality  inconsistency  with  this  form  of  weak  realizability  that  led  us  to  propose 
an  alternative  form  which  allows  model  predictions  to  come  ubitrarily  close  to  one 
and  two  component  states.  Interestingly  enough,  the  Fu,  Launder  and  Tselepidalds 
(FLT)  model  satisfies  this  alternative  form  of  weak  realizability  and  thus,  unlike  the 
Shih-Lumley  model,  it  guarantees  positive  component  energies  in  general  homogeneous 
turbulent  flows. 

(4)  In  an  attempt  to  satisfy  the  strong  form  of  realizability,  both  the  Shih-Lumley  and  FLT 
models  introduce  higher  degree  nonlinearities  in  the  modeling  of  the  rapid  pressure- 
strain  correlation.  It  was  demonstrated  in  computations  of  homogeneous  shear  flow 
that  these  higher-order  terms  cause  the  models  to  become  highly  ill-behaved  near  one 
or  two  component  states  of  turbulence  where  they  generate  large  amplitude  oscillations 
that  are  unphysical.  Older  models  such  as  the  Launder,  Reece  and  Rodi  and  IP  models 
that  are  linear  do  not  have  this  problem.  It  must  be  remembered  that  the  rapid 
pressure-strain  correlation  is,  by  its  definition,  linear  in  the  energy  spectrum  tensor  - 
a  property  that  such  nonlinear  models  do  not  possess. 
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Finally,  some  comments  are  warranted  concerning  the  important  implications  that  these 
results  have  for  turbxilence  modeling.  The  notion  that  realizability  constraints  can  be  used 
to  calibrate  a  second-order  closure  model  appears  now  to  be  fundamentally  unsound.  The 
process  of  ensuring  realizability  is  decidedly  non-unique.  There  is  no  apparent  way  to  satisfy 
the  strong  form  of  realizability  in  the  existing  hierarchy  of  second-order  closures  and  there 
are  an  infinity  of  ways  in  which  the  weak  form  of  realizability  can  be  implemented. 

This  brings  us  to  the  more  fundamental  question:  Why  is  there  such  an  obsession  with 
realizability  constraints  in  second-order  closure  modeling?  It  is  a  weU  established  result 
that  the  standard  form  of  the  modeled  dissipation  rate  equation  (12)  guarantees  limited 
realizability  -  namely,  non-negative  values  for  the  turbulent  kinetic  energy  and  dissipation 
rate  in  homogeneous  turbulence  (see  Speziale  1990).  Furthermore,  so  long  as  the  fixed  points 
of  the  model  are  realizable  and  sufficiently  strong  attractors,  at  worst  there  will  be  a  short 
transient  where  one  or  two  components  of  the  turbulent  kinetic  energy  become  negative. 
This,  in  general,  does  not  have  to  be  computationally  fatal,  particularly  if,  for  inhomogeneous 
flows,  the  turbulent  diffusion  terms  are  modeled  with  isotropic  gradient  transport  models 
such  as  those  introduced  by  Mellor  and  Herring  (1973)  (numerical  instabilities  generally  arise 
from  negative  diffusivities).  Typically  the  short  transients  where  realizability  violations  occur 
are  in  turbulent  flows  that  are  so  far  firom  equilibrium  that  simple  one-point  closures  cannot 
be  expected  to  apply  in  the  first  place.  In  an  effort  to  avoid  short  lived  regimes  where 
the  solution  can  become  unrealizable,  models  developed  based  on  realizability  constraints 
have  been  complicated  substantially,  leading  to  highly  ill-behaved  solutions  that  have  no 
greater  predictive  capabilities  than  those  of  the  simpler  models  near  the  limits  of  realizable 
turbulence.  In  fact,  the  higher-degree  nonlinearities  introduced  to  satisfy  realizability  often 
render  a  model  ill-behaved  to  the  point  where  the  solution  is  contaminated  by  large  amplitude 
oscillations  for  long  durations  of  time.  This  is  extremely  unwise.  It  would  be  far  preferable 
to  introduce  mathematical  devices  to  avoid  computationally  dangerous  unrealizable  behavior 
in  turbulent  flows  that  are  far  from  equilibrium  without  compromising  the  near-equilibrium 
predictions  of  a  model.  A  simple  mathematical  means  for  making  any  existing  second-order 
closure  model  realizable  in  this  fashion,  based  on  a  stochastic  analysis,  is  the  subject  of  a 
companion  paper  (Durbin  and  Speziale  1993). 
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APPENDIX  A 


Relative  to  a  coordinate  system  with  unit  vectors  Aq  that  is  undergoing  a  time- dependent 
rotation,  any  symmetric  tensor  T  takes  the  form: 

T  =  TotfiXaXff.  (-^1) 

Hence,  the  first  time  derivative  is  pven  by 

T  =  TafiXaXfi  +  TaffXaXff  -h  Tq^AqA^  (-^2) 

since  \q  =  Aa(t).  It  is  wed  known  that  for  any  unit  vector  Aq  rotating  with  angular  velocity 
n  we  have  (cf.  Goldstein  1980) 

=  nxXa 


Hence, 


or  equivalently 


—  ^Aya^yXs. 


T  =  {Ta0  +  eayS^yTs0  -j-  e0^s^yTsa)XaX0 


(T)a;8  =  Ta0  +  CaygHyTsff  -f 


The  second  time  derivative  is  given  by 


T  =  TaffXaXff  -|-  2Ta0XaXp  -f-  2Ta0XaX0  -f-  2Ta0XaXfi  -f  TaffXaXff  -f  TaffXaXff.  (■<46) 


Since, 


Aa  =  DxAa-l-DxAo 


we  have 


=  i7xAa -J- Dx(DxAq) 

=  eSya^yXg  -}-  ila^gXg  ~  fi*Aa 
=  {Ta0  +  2ecr7«n^!ZV/3  -f  2e^jsiljTsa 
'^^^ayS^0i»T^-y^i/1’sir  "I"  ^a’f6^s'I's0 


^0y6^yTsa  +  ^<S^sTs0  +  ^0^sTso 
-2n^Ta0)X^X0 


where 
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Capelli’s  identity,  given  by  the  determinant 


^-jr^  ^Tr<r 

Sx^  Sxu  8x„ 


can  be  applied  to  the  fourth  term  on  the  right-hand-side  of  (A8).  After  simplifying, 
at  the  final  result: 

(^)a/9  ~  'fa$  "I"  2eoi-Yjfl-jT]5(J  -f-  "I"  ^orysiiyTgff 

-f  -|-  ZQa^rT,T0  +  ZQffClffTaa 

-4Q^Tc0  -  2T„„Qa^li  +  -  n^n.T,<r)^a^. 


(A9) 

arrive 


(AlO) 
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APPENDIX  B 


The  detailed  form  of  the  turbulence  models  considered  in  the  paper  are  as  follows: 
Shih  &  Lumley  Model 


Dij  =  —Ciebij  4-  —K^ij  +  12a5A'  {hn^jk  +  bjkSit 
2  _  \  4 

—  ^buSuSijj  +  -(2  —  7as)K{bik<ijjk  +  bjti^ik) 

+  jK{bnblmSjm  +  bjlblmS  im  ^bik^kibij 

5 


(Bl) 


^bklSklbij)  4"  ^K{bilblfnU jm  "I"  bjlblfnt*fitn) 

C,i  =  1.20,  Cc2  =  l  +  0.49 exp(-2.83//^)[l  -  0.33/n(l  -  55//)] 

0 

2  dxi 

Cl  =  2  +  ~exp(-7.77/y^){72/y^  4-  80.1  /n[l  +  62.4(-//  4-  2.3///)]} 

B  =  1  4-  9//  4-  271 1 1 

II=-\bijbij,  III  =  ^bijbjkbki 


Ret 


4^ 
9  ve 


1 


/h,  Launder  &  Tselepidakis  Model 


no 


—Cicbij  4- 


(B2) 

(B3) 

m 

(B5) 

(B6) 

(B7) 

(B8) 


4-  1.2JK’  (bik^jk  +  bjk^ik  -  ^bkiEuSij^ 

26  4  — 

+TzK{bikUjk  4-  bjklDik)  +  ■zK{bikbkiSji 
15  5 

+bjkbkiSii  —  2bikSkibij  —  2bkiSkibij) 


(B9) 


4  14  _ 

4--A  {bikbki(tfji  4-  bjkbkiu/ii)  — —K  [BII{bik(^jk 

+bjktOik)  +  I2{bikbia0flmbmj  +  bjkbuH^ilmbmi)] 
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Cel  =  1.44,  Ce2  =  1.90 

Cl  =  -  -  2Vf  +  2,  Cj  =  uaiiVf 


(J510) 

(511) 


IP  Model 


where 


n,j  —  Cichij 

Cel  =  1.44,  Ce2  =  1.90 


Cl  =  3.6,  C2  =  0.6 

_  dvj  dvi 

r,j  —  Tile  Q  Tjic  ^ 

oxic  axie 


(fll2) 

(B13) 

(B14) 

(B15) 


■P=-Ti 


•’8*/ 


(B16) 
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Figure  1.  Schematic  of  possible  time  evolutions  of  the  invariant  function  F  passing 
the  realizability  boundary  F  =  0. 


(a)  SHIH-LVMLEY  MODEL 


-0.04  0.00  0.04  0.08 


III 

Figure  4.  Computed  phase  space  trajectories  in  homogeneous  shear  flow  for  the  initial 
conditions  (6ii)o  =  —0.32,  (622)0  =  —0.32,  SKo/eq  =  50:  (a)  Shih-Lumley  Model  and  (b) 
the  IP  Model  of  Launder  and  co-workers. 
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Figure  4.  Computed  phase  space  trajectories  in  homogeneous  shear  flow  for  the  initial 
conditions  (6ii)o  =  —0.32,  (622)0  =  —0.32,  SKq/co  =  50:  (a)  Shih-Lumley  Model  and  (b) 
the  IP  Model  of  Launder  and  co-workers. 
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Figure  5.  Model  predictions  for  the  time  evolution  of  the  Reynolds  stress  amsotropies  in 
homogeneous  shear  flow  corresponding  to  the  initial  conditions  (feii)o  =  —0.32,  (622)0  = 
—0.32,  SKofeo  =  50:  (a)  F  and  (b)  612. 
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Figure  5.  Model  predictions  for  the  time  evolution  of  the  Rejmolds  stress  anisotropies  in 
homogeneous  shear  flow  corresponding  to  the  initial  conditions  (611)0  =  —0.32,  (622)0  = 
-0.32,  SKoleo  =  50:  (a)  F  and  (b)  612. 
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Figure  6.  Computed  results  for  homogeneous  shear  flow  corresponding  to  the  initial  con¬ 
ditions  (hii)o  =  —0.3,  (622)0  =  —0.3,  SKofeo  =  10:  (a)  phase  space  trajectories  of  the 
Shih-Lumley  Model,  (b)  time  evolution  of  F  and  (c)  time  evolution  of  612. 
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Figure  6.  Computed  results  for  homogeneous  shear  flow  corresponding  to  the  initial  con¬ 
ditions  (611)0  =  -0.3,  (622)0  =  -0.3,  SKofeo  =  10:  (a)  phase  space  trajectories  of  the 
Shih-Lumley  Model,  (b)  time  evolution  of  F  and  (c)  time  evolution  of  612. 
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Figure  6.  Computed  results  for  homogeneous  shear  flow  corresponding  to  the  initial  con- 
ditions  (611)0  =  —0.3,  (622)0  =  —0.3,  SKofeo  =  10:  (a)  phase  space  trajectories  of  the 
Shih-Lumley  Model,  (b)  time  evolution  of  F  and  (c)  time  evolution  of  612. 
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Figure  7.  Computed  phase  space  trajectories  for  homogeneous  shear  flow  corresponding  to 
the  initial  conditions  (611)0  =  —0.32,  (622)0  =  —0.12,  SKo/eo  =  50:  (a)  the  Shih-Lumley 
Model  and  (b)  the  IP  Model  of  Launder  and  co-workers. 
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Figure  7.  Computed  phase  space  trajectories  for  homogeneous  shear  flow  corresponding  to 
the  initial  conditions  (6ii)o  =  —0.32,  (622)0  =  —0.12,  SKoIeq  =  50:  (a)  the  Shih-Lumley 
Model  and  (b)  the  IP  Model  of  Launder  and  co-workers. 


Figure  8.  Computed  phase  space  trajectories  for  homogeneous  shear  flow  corresponding  to 
the  initial  conditions  (611)0  =  -0.24,  (622)0  =  0.17,  (612)0  =  0.2,  5ifo/eo  =  50:  (a)  the  IP 
Model  and  (b)  the  FLT  Model  of  Launder  and  co-workers. 


Figure  8.  Computed  phase  space  trajectories  for  homogeneous  shear  flow  corresponding  to 
the  initial  conditions  (611)0  =  -0.24,  (622)0  =  0.17,  (612)0  =  0.2,  SKo/eo  =  50:  (a)  the  IP 
Model  and  (b)  the  FLT  Model  of  Launder  and  co-workers. 
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Figure  9.  Model  predictions  for  the  time  evolution  of  the  Reynolds  stress  anisotropies  m 
homogeneous  shear  flow  corresponding  to  the  initial  conditions  (611)0  —  0.24,  (622)0  0.17, 

(612)0  =  0.2,  and  SKofso  =  50:  (a)  F  and  (b)  612. 
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Figure  9.  Model  predictions  for  the  time  evolution  of  the  Reynolds  stress  amsotropies  m 
homogeneous  shear  flow  corresponding  to  the  initial  conditions  (hii)o  =  —0.24,  (622)0  =  0.17, 
(612)0  =  0.2,  and  SKofco  =  50:  (a)  F  and  (b)  612. 
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F'igure  10.  Computed  phcise  space  trajectories  of  the  FLT  Model  in  homogeneous  shear  now 
:orresponding  to  the  initial  conditions  (6ii)o  =  —0.32,  (622)0  —  —0.32,  SKq/eq  =  50. 


IP  MODEL 


(o) 


0  25  50  75  100 

t* 


Figure  11.  Computed  model  predictions  for  the  time  evolution  of  the  ratio  of  production  to 
dissipation  (P/s)  in  homogeneous  shear  flow  corresponding  to  the  initial  conditions  (hn)o  = 
—0.32,  (622)0  =  —0.32,  SKo/eo  =  50:  (a)  IP  Model,  (b)  Shih-Lumley  Model  and  (c)  the  FLT 
Model. 
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Figure  11.  Computed  model  predict! 
dissipation  {Vje)  in  homogeneous  sh 
-0.32,  (622)0  =  -0.32,  SKolto  =  50: 
Model. 
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Figure  12.  Time  evolution  of  the  normal  Reynolds  stress  r22/(r22)o  in  homogeneous  shear 
flow:  Comparison  of  computed  results  for  the  FLT  Model  (corresponding  to  (&ii)o  =  (^22)0  = 
—0.0833,  SKoleo  =  50)  with  the  Rapid  Distortion  Theory  (RDT)  solution. 
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